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Abstract 

The rate of escape of an ideal bead-spring polymer in a symmetric double-well potential 
is calculated using transition state theory (TST) and the results compared with direct 
dynamical simulations. The minimum energy path of the transitions becomes flat and the 
dynamics diffusive for long polymers making the Kramers-Langer estimate poor. However, 
TST with dynamical corrections based on short time trajectories started at the transition 
state gives rate constant estimates that agree within a factor of two with the molecular 
dynamics simulations over a wide range of bead coupling constants and polymer lengths. 
The computational effort required by the TST approach does not depend on the escape 
rate and is much smaller than that required by molecular dynamics simulations. 
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1 Introduction 


The movement of a polymer from one metastable state to another (meta)stable state by 
crossing a free energy barrier through thermal activation is referred to as the polymer escape 
problem. The free energy barrier separating the states is most often of entropic origin due 
to geometric confinement. Relevant examples include polymer translocation mm, where a 
polymer is crossing a membrane through a pore [3], or narrow /xm-scale channels with traps 
[1]. Recent experiments by Liu et al. involve the escape of DNA molecules from an entropic 
cage [5]. Similar translocation and escape processes are common in cell biology and have 
possible bioengineering applications, such as DNA sequencing [B] and biopolymer filtration 

m- 

The polymer escape problem is a generic description of systems such as the ones described 
above. In the polymer escape problem there is a static barrier represented by an external 
potential energy function where the two minima are equal or the final state has lower energy 
than the initial state. Therefore, the escape problem is a combination of the classical thermal 
activation problem, where the free energy barrier is largely energetic urn, and the polymer 
translocation problem, where the free energy barrier is typically of entropic origin. 

The problem has been studied using Kramers’ theory for polymers by Park and Sung, 
who evaluated the escape rate using lattice statistics of a discrete ideal polymer model [TO] . 
Sebastian proposed a kink diffusion mechanism m for long chains when one end of the 
polymer has moved over the barrier, while the other end is still in the initial state energy 
well. The kink corresponds to the beads in the region of the energy barrier and it moves 
along the chain as the polymer moves from one potential energy well to another. In Ref. 
m Sebastian and Paul described the polymer escape of long chains as a two-step process. 
In the first step, the polymer is thermally activated to bring one end over energy barrier 
and into the final state well, and the second step involves diffusive motion of the kink as 
intermediate beads move in and out of the barrier region. A rate theory approach similar 
to hanger ’s m multidimensional extension of Kramers’ theory was proposed for activation 
m- More recently, Sebastian and Debnath studied the thermal activation mechanism for 
short chains m and simulated kink diffusion for one and three dimensional systems |15] . Lee 
and Sung studied the polymer escape in a symmetric external potential well and proposed a 
rate theory approach to predict the rate for linear m and for star polymers m- They also 
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found that for linear polymers the stretched kink solution is the dominant escape mechanism 
for chains longer than a certain crossover length Nc- Below Nc the polymer crosses the 
barrier in a coiled form while polymers that are longer than Nc are stretched during the 
transition, analogous to instantons in quantum mechanical tunnelling of one particle. Paul 
|18] has studied polymer escape of star polymers in a system mimicking experiments carried 
out by Han et al. [1]. 

The studies mentioned above only considered the escape of flexible, ideal polymers with 
zero bond length at zero temperature. The semiflexible case has been studied by Kraikivski 
m- Self-avoiding (SA) polymer models have been studied numerically and compared with 
the flexible and semiflexible ideal chain models in Refs. [20] and mi SA polymers show 
qualitatively similar behaviour as ideal polymers, but unlike the monotonically decreasing 
escape rate of ideal polymers, the escape rate of the SA polymers exhibits a minimum for 
intermediate length beyond which the escape rate increases. 

In this article, we study polymer escape dynamics and estimate the escape rate using 
harmonic approximation of the transition state theory followed by dynamical corrections 
(DC) for a discrete, one-dimensional harmonic ideal polymer model in a symmetric double¬ 
well external potential. We compare the results with molecular dynamics simulations using 
both Brownian dynamics (BD) and Langevin dynamics (LD). We also compare with results 
obtained from hanger’s rate estimate [mE]. Using the Nudged Elastic Band (NEB) method, 
we find the minimum energy paths (MEP) to identify the relevant saddle points on the energy 
surface as a function of chain length N and polymer spring constant K. The harmonic modes 
at the saddle point are evaluated and their eigenvalues and eigenvectors used to analyse the 
escape dynamics of the polymer. 

The system studied is described in Sec. [21 the methods for calculating the rate in Sec. 
E] and the numerical methods in Sec. [H The results are presented in Sec. EJ The article 
concludes with a summary and discussion in Sec. [6l 

2 System Description and Problem Definition 

We consider a polymer represented by N beads in one dimension interacting with their neigh¬ 
bours through a spring. The configuration of the polymer is described by the set of coordinates 
r := {f’n}n=i with the centre of mass Yln=i '‘"n- The equation of motion for the nth bead 
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at time t is given by Langevin dynamics (LD) as 


mrn{t) + (t) + '^n[V irn{t)) + U] = y^2-fkBT^n{t), (1) 

where m is the mass of an individual bead, 7 the friction coefficient, V{rn) the external 
potential, U the interaction potential between beads, ksT the thermal energy, and Cn{t) a 
Gaussian random force satisfying {Cn{t)) = 0, and {Cn{t)Cm{t')) = 6{t — t')5n,m- In the limit 
of strong coupling with the heat bath, i.e. in the high friction limit, the dynamics become 
overdamped and are described by the Brownian equation of motion (BD) 

7r„(t) + Vn[V irn{t)) + U] = y/2qk^^n{t)- (2) 


The interaction potential between the beads is given by a harmonic potential U = 
- Tn+if giving the force 

-VnU =-K{rn-i + rn+i-2rn), (3) 


acting on the nth bead. The external potential V{x), shown in Fig. [H is a quartic double 
well 

where ±ao gives the location of the minima, the energy has a maximum at x = 0 , and 
is the curvature of the energy barrier. The same polymer model and external potential with 
BD were used in Ref. M The energy of a polymer configuration r is given by 


N 


N-1 


^(r) = “ rn+if. 


(5) 


n=l 


n=l 


In the one-dimensional Rouse model, a continuum model of an ideal chain, the spring 
constant is iF = ksT/l^ which gives the polymer an average bond length /q = y/ksT/K 
and an average squared end-to-end distance of (Rgg) = NIq in free space. The parameter 
Roc = \/ (Ree) can be used to compare the size of the polymer to the width of the potential 
well. [ 22 ] 

At time t = 0, all the beads of the polymer are in the left potential well shown in Fig. [H 
so Aq < 0. The polymer is thermally equilibrated with the heat bath and its configuration 
consistent with Boltzmann distribution F’i(r) = P(r|Ao < 0) oc exp(— d>(r)/A;BT). With time 
evolution according to either Eq. ([T]) or Eq. ([2|), the polymer eventually escapes to the final 
state E, where Xq > 0. The rate of such escape events is TZi-^p and the time in between 
escape events is tcross ~ average. 
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3 Calculation of the Rate 


3.1 Molecular Dynamics Simulations 


The LD Eq. ([T]) or the BD Eq. ([2|) are integrated numerically to generate trajectories 
representing transitions from I to E. Uncorrelated samples of the Boltzmann distribution 
Pi(r) for the polymer in state I were generated by simulating the chain confined to the left 
potential well for time intervals twice as long as the correlation time between samples. Erom 
these trajectories, the escape probability was determined as 

Pesc(i) = (l/iVtraj) (6) 


2 = 1 

where A^traj is the number of simulated trajectories, 9{t — U) the Heaviside step function and 
ti the time of the ith escape event. An escape event is considered to have occured when the 
centre of mass of the system Xq has reached Xq > ao/2. The dynamical escape rate is given 
by 




dPescjt) 

dt 


(7) 


where the derivative is computed by fitting the curve Pesc = P-MDt + h, where b is a constant, 
over the time interval where Pesc(i) is close to being a linear function of t. The average time 
in between escape events is tcross = (U) = (1/Atraj) 


3.2 Transition State Theory 

Transition state theory (TST) assumes that the initial distribution has the Boltzmann form 
Pi(r) oc exp(—<h(r)//cBT) and that the energy barrier is high enough for the time between 
escape events to be longer than the relaxation time. A transition state is defined in between 
the initial and final states that should represent a bottleneck for the transition. The key as¬ 
sumption of TST is the no-recrossing approximation where a trajectory crossing the transition 
state with velocity pointing away from the initial state is assumed to end up in the final state 
without recrossings. Implicit in this approximation is an assumption that the time spent in 
the vicinity of the transition state is short compared to the time spent in the final state. The 
TST rate estimate is obtained by multiplying the probability that the system reaches the TS 
with the flux through the TS 

T^tst = \/kBT/{2Trfi±)Z'Ts/Zi, (8) 
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where fj,± is the reduced mass of the system in the direction perpendicular to the TS surface, 
and Z'j's and Zi are the configuration integrals of the TS and I, respectively. [8] 

Instead of computing the full configuration integrals Zts and Zi, a harmonic approxima¬ 
tion to TST (HTST) can be used [23]. If the escape trajectories are likely to go through 
the vicinity of a first order saddle point on the energy surface in between the energy minima 
corresponding to the I and F states, the configuration integrals can be evaluated analytically 
by making a second order expansion around the saddle point and the initial state minimum. 
A convenient choice for a reaction coordinate is the minimum energy path (MEP) and the 
relevant saddle point is the point of maximum energy along the MEP. In HTST, the TS is 
chosen to be a hyperplane going through the saddle point with its normal pointing along the 
MEP. The MEP starts from rj = —oqI, where 1 = [1,1,... , 1] is a vector of length N, and 
ends at rp = oq! and the location of the saddle point is dennoted by r^. After expanding the 
total energy of the system up to second order around the initial state minimum tq and the 
saddle point rj, the configuration integrals Zts and Zi can be evaluated analytically to give 
the HTST estimate of the transition rate 


T^htst 


1 

2TTy/JIZ 



A<S>/kBT 


(9) 


Here A? and A| are the eigenvalues of Hessian matrix which has elements 


{Hp)ij — 


(92$(r) 


dvidrj 


( 10 ) 


evaluated at the initial state minimum and at the saddle point. The product in the de¬ 
nominator of the HTST rate expression omits the negative eigenvalue corresponding to the 
direction along the MEP. The elements of the Hessian can be computed using finite differences 
of the forces acting on the beads. Diagonalisation of the Hessian matrix of Eq. m gives an 
eigenvalue spectrum diag(Aj) = S“^HS and a set of eigenvectors, S = [si ... sjv], which are 
referred to as the harmonic modes. 

The configuration of the polymer at the saddle point depends on the chain length in such 
a way that below a crossover length Nc the polymer is collapsed to one point r| = 0. For 
longer polymers with N > Nc, it becomes energetically more favourable for the chain to be 
stretched along the MEP. This transition from collapsed to stretched polymer is analogous to 
the onset of quantum mechanical tunneling and the appearance of an instanton in quantum 
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rate theories based on Feynman path integrals. The relation between the barrier curvature, 
spring constant and critical chain length can be derived analytically from the continuum 
model of an ideal chain m- 

I Ktt^ 

At Nc, the eigenvalue of the smallest positive mode at the saddle point, A^, approaches 
zero. This leads to a divergence in the HTST rate estimate which can be reduced by introduc¬ 
ing anharmonic corrections (AHC) for this mode. Lee and Sung have derived an expression 
for the correction factor as 

»(“) = /” (12) 

where a = aIoq y/N/{kBT)/uj [16]. Thus, the corrected HTST rate becomes 

^HTST+AHC = 5(«)'^HTST- (13) 

3.3 Dynamical Corrections 

The assumption of no recrossings of the TS in TST is often not satisfied for energy barriers 
that are broad. The effect of the recrossings can, however, be estimated using calculations of 
classical trajectories started at the TS. Initial configurations of the trajectories are sampled 
within the TS hyperplane, such that net force is zero Fq = F — Fi/N = 0 which keeps 
the system within the plane. This provides dynamical corrections (DC) to the TST rate 
estimate. The correction factor, n, is |24] : 

j;sgn(u,)0(Afi-i), (14) 

where sgn(ui) is the sign of the initial net velocity Vi = ^ trajectory i assigned 

from the Maxwell-Boltzmann distribution P{v) (x exp[—u^/(2/c5r)], and is unity if 

the system resides in F at the end of the trajectory and zero otherwise. Thus, each trajectory 
that ends in F, contributes positively to k if it starts with a velocity pointing towards F, but 
contributes negatively if the initial velocity points towards I. Trajectories ending in I do not 
contribute to n. The corrected rate estimate is 

'^HTST+AHC+DC = «^^HTST+AHCi (15) 
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where T^htst+ahc is given by Eq. (fT^ . While the direct dynamical calculations of transitions 
can be impossibly long for systems with low transition rates, the calculation of DC requires 
only short trajectories and can be carried out readily with little computational effort. 


3.4 Kramers’ and Danger’s Rate Theory 


hanger’s rate theory for barrier crossing is a multidimensional generalisation of Kramers’ 
canonical rate theory m- hanger’s expression for the rate resembles the HTST rate expres¬ 
sion, although it is derived from different assumptions. The main difference between hangers’ 
and HTST rate estimates is that hanger’s rate has an additional pre-exponential factor that 
provides an approximate estimate of the effect of TS recrossings, hanger’s rate estimate is: 


7^L = 




A^/kgT 


(16) 


2>r7 \| n£, A* 

hanger’s approach was used by Sebastian et al. and hee and Sung They calculated 

the unstable mode and partition functions analytically for the Rouse model, a continuum 
model of an ideal chain, hanger’s rate estimate diverges around Nc similar to the HTST rate 
estimate and this divergence can be removed using the correction factor of Eq. m- 


4 Numerical Methods 

The escape rate was obtained by numerical molecular dynamics simulations using both 
hangevin dynamics (ED) of Eq. ([1]) and Brownian dynamics (BD) of Eq. ([2|). BBK in¬ 
tegration [2S] was used for hD with a time step of At = 0.005, and forward Euler integration 
for BD with time step ranging between 0.005 and 0.01. Chains with different combinations 
of N and K in ranges N £ {1,... , 120} and K £ {5,... , 60} were simulated with parameters 
7 = 1.0, m = 1.0 and ksT = 1.0. The parameters for the external potential of Eq. (BD were 

= 1.5 and Uq = 1.5. The parameters are the same as those used in Ref. [161 If we choose the 
units of length, mass and energy to be Iq = 1.02 nm, mo = 1870 amu, corresponding to double 
stranded DNA, and ksT at T = 300 K, the unit of time becomes to = -\/ molQ/ksT = 27.9 
ps. 

In order to find the relevant saddle point to apply HTST of Eq. Q the Nudged Elastic 
Band (NEB) method [26l was used to find the MEP of the escape transition. In the NEB 
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method, a set of replicas of the system, referred to as images, {rp}^=i are placed along a 
path between the initial state minimum ri = —aol and the final state minimum rp = aol- 
The images represent a discretisation of the path and to control the distribution of these 
discretisation points, the images are connected with harmonic springs with spring constant 
A:neb • An estimate of the tangent m to the path is used to project out the parallel component 
of the true force acting on the beads and the perpendicular component of the spring force 
(the ’nudging’). The images are then displaced iteratively so as to zero the net force acting 
on the beads until the images lie along the MEP. 

At the end of the NEB calculation the images give a discrete representation of the MEP, 
but no image will be exactly at the saddle point. An accurate estimate of the saddle point 
was found here by identifying the image with the highest energy and then minimising the 
force on this image with the Newton-Raphson method. The method consists of iterations 
j,(n+i) _ j.(n) _ where F®y®(r) is the force of the real system (no NEB spring 

forces) and H~^ is an inverse of the Hessian matrix (Jacobian matrix of the force) with 
elements given by Eq. m- Within a few iterations, the method converges to a point rj 
where F(rj) = 0. If the initial point is close enough to a saddle point, then the Newton- 
Raphson method will most likely converge onto the saddle point rather than other points 
where the force is zero. The activation energy is then given by the energy difference between 
the initial state and the saddle point as A<I> = <h(r|) — <k(ro). The Hessian matrices Hq and 
H| of Eq. (HOD can be evaluated using the finite difference method. Alternatively, the exact 
saddle point can be found with the climbing image method [28]. 

The spring constants used in the NEB calculations was A:neb = 8.2 x 10“® for K = 10, 
and ^neb = 8.2 for K = 60. The number of images, P, was typically chosen to be between 9 
and 19, but for the larger JV sometimes up to P = JV/2. 

5 Results 

5.1 Minimum Energy Paths and Saddle Points 

The energy along minimum energy paths for polymer escape are shown in Fig. [2ja. The 
energy barrier is found to increase linearly with N up to Nc (the integer value of the crossover 
length obtained from MEP calculations is denoted by Nc while the estimate obtained from 
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the continuous ideal chain, Eq. dm), is denoted by Nc) where it reaches a constant value. 
Below Nc the saddle point configuration corresponds to a coiled up polymer where all the 
beads lie on top of the maximum of the external potential at x = 0, and the height of the 
energy barrier is A'h = NAV. Beyond Nc, the stretched configuration is energetically more 
favourable. The values of Nc obtained from MEPs are Nc = 8 for K = 10 and Nc = 20 for 
K = 10, which agree well with theoretical predictions of Eq. (jll|) : Nc = 8.11 for K = 10 
and Nc = 19.86 for K = 60. 

The eigenvalue of the smallest positive mode at the saddle point, A^, decreases with N 
until Nc, where it has a minimum of nearly zero. Beyond Nc it increases, as shown in 
Fig. [3la. The mode corresponding to this eigenvalue is responsible for creating the stretched 
configuration; instead of the eigenvalue aI becoming negative and the mode unstable, the 
polymer stretches along the energy barrier of the external potential, extending the springs 
between beads near the barrier top but at the same time lowering the number of beads at the 
barrier top and thereby reducing the overall activation energy below NAV. The stretched 
configuration for the chain of length A = 56 is shown in Fig. [3lb along with the modes 
corresponding the eigenvalues A| and aI . 

The mode with the negative eigenvalue A| corresponds to the polymer moving towards 
the final state as shown in Fig. Ob. Below Nc, the negative eigenvalue A| corresponds to the 
curvature of the external potential At Nc, the negative eigenvalue A| starts increasing 
and approaches zero with increasing N, because the energy along the MEP flattens out. Fig. 
Ob shows the energy along the MEP for various values of N illustrating how the curvature 
at the saddle point decreases with increasing chain length. 

5.2 Brownian and Langevin Dynamics 

The molecular dynamics simulations were performed using both Langevin dynamics of Eq. 
(dl) and Brownian dynamics of Eq. ©• An escape event was considered to have occured when 
the centre of mass of the system Xq had reached Xq = ao/2 ~ 0.61. As can be seen in Fig. 
Ob, the barrier is flat in this region for the longer chains. We repeated simulations for chains 
N G {56, 64, 80} with K = 10 and the crossing condition Xq = 1.0, and found the difference 
in the calculated rates to be negligible. 

The escape rate obtained from these simulations using Eq. Q is shown in Fig. 0 as a 
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function of N for chains with K = 10 and K = 60 and as a function of the spring constant 
K for chain of length = 64 in Fig. [5j The maximum chain length for for which LD 
simulations were carried out was limited by computational resources. The escape trajectories 
were started from the equilibrated initial distribution T’i(r) oc exp{—^{r)/kBT) with the 
constraint Xq < 0 for both LD and MD. Lee and Sung used the initial condition r = — ool 
m instead of trajectories started from a Boltzmann distribution. We compared these two 
initial conditions and found no difference in the escape rate. 

Both escape rates show a monotonically decreasing rate in both N and K, which is 
consistent with previously published work [Ml 123 E]. Our estimates of the escape rate 
obtained from BD agree well with the results reported by of Lee and Sung |16] . The rate 
obtained from BD trajectories is slightly higher than the rate obtained from LD trajectories. 
We note that the BD and LD rates are not the same because of the lack of inertial effects 
in BD. Nevertheless, with the current parameters the LD simulation results are very close to 
the overdamped BD case. 

5.3 Transition State Theory with Dynamical Corrections 

The escape rate 7^htst+ahc obtained by harmonic transition state theory (HTST) with 
anharmonic corrections (ABC) of Eq. (fMjl - and 7^htst+ahc+dc with dynamical corrections 
(DC) of Eq. ()15p . are shown in Fig. [4] as a function of N for chains with spring constants 
K = 10 and K = 60. In Fig. [5l T^htst+ahc and 7^htst+ahc+dc are shown as a function of 
the spring constant K for a chain of length N = 64. The peak at Nc arising from Eq. ([3) is 
removed from the rate in Fig. fusing the anharmonic correction of Eqs. (1121) and (I13p . The 
HTST rate of Eq. Q without AHC is compared with the corrected rate in Fig. [6j The peak 
arises in the HTST rate expression because the eigenvalue of the smallest positive mode at 
the saddle point approaches zero as shown in Fig. [3la. For the continuum chain model 
without anharmonic corrections this peak would diverge to infinity as Ag —)• 0 at Nc- 

Beyond Nq the saddle point bifurcates into two saddle points where the chain is stretched 
over the barrier of the external potential as shown in Fig. [3jb. There are two saddle points 
corresponding to the stretched configuration since either one of the two ends can be displaced 
towards the final state. Since the chain can escape through either of these symmetric saddle 
points, the total rate is twice the HTST rate given by Eq. ([9]) when N > Ne¬ 
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The rate estimate obtained by HTST with AHC saturates to a constant value quickly 
beyond Nc as shown in Fig. SI This is because the energy rise along the MEP as well as the 
positive eigenvalues reach constant values. However, the rate obtained from the molecular 
dynamics simulations continues to decrease with N. The energy along the MEP becomes 
flat as illustrated in Pig. [21b leading to increased number of recrossings of the TS. When DC 
obtained from short time classical trajectories started at the TS are included, the rate estimate 
of Eq. (1151) . close agreement is obtained with the rate obtained from the direct molecular 
dynamics simulations, the latter being typically a factor of two larger for the longest polymers 
studied. 

5.4 Danger’s Theory 

The escape rate given by Langer’s theory of Eq. (I16p is shown in Pig. |4] as a function of N 
and in Pig. [2]as a function of K. Eor short polymers, below Nc-, Langer’s rate estimate agrees 
well with the rate obtained from the molecular dynamics simulations but it underestimates 
the rate by orders of magnitude for the long chains. The reason for this is the flattening of 
the energy along the MEP which causes the negative eigenvalue, A|, corresponding to the 
unstable mode at the saddle point, to approach zero. The Langer-Kramers estimate of the 
recrossing correction which results in the appearance of A| in the prefactor in Eq. (1161) then 
gives a large overcorrection and much too small rate. 

6 Summary and Discussion 

The results presented here on the polymer escape problem for the one-dimensional bistable 
external potential of Lee and Sung m show that a TST approach based on calculations 
of MEPs to identify saddle points corresponding to the escape transitions followed by a 
HTST rate estimate with anharmonic and dynamical corrections gives good agreement with 
results obtained with direct MD simulations for a wide range in polymer length and spring 
interaction between beads. These TST calculations require little computational effort and can 
be extended easily to three-dimensional systems with complex interaction potentials, while 
the MD simulations are limited to short polymers and simple interactions. The computational 
cost of finding saddle points and computing recrossing trajectories does not depend on the 


12 


activation energy for the escape. The flat energy profile along the MEP in the diffusion 
regime makes the dynamical correction trajectories somewhat longer but they are still orders 
of magnitude shorter than direct MD simulations of the escape transitions. The harmonic 
mode analysis, furthermore, gives valuable insight into the transition mechanism. The two 
lowest harmonic modes at the saddle point, shown in Fig. [3l explain the relevant escape 
dynamics of the chain. The lowest, unstable mode corresponds to flux along the MEP towards 
the hnal state while the mode corresponding to the second lowest eigenvalue, causes the chain 
to become stretched for N > Nc, and the saddle point to bifurcate. The eigenvalue of the 
unstable mode approaches zero as N increases further leading to a flat energy profile along 
the MEP, as was found by Lee and Sung m- This makes the escape process diffusive along 
the MEP and the system recrosses the transition state multiple times, calling for explicit 
dynamical corrections to HTST using short time trajectories started at the TS. The ca. 
factor of two underestimate of the rate for long chains is likely due to an underestimate of the 
entropy of the transition state with respect to the entropy of the initial state in the harmonic 
approximation. If a full TST calculation were carrier out, followed by a calculation of the 
dynamical correction, the rate estimate would be exact and agree with the direct dynamical 
simulations. 

The rate estimate obtained from the Kramers-Langer approach which has previously been 
used extensively [HI [HI [H] turns out to give a poor estimate for the longer chains. The 
reason is the flat energy profile along the minimum energy path in the diffusion regime which 
causes the eigenvalue corresponding to the unstable mode at the saddle point to approach 
zero. As a result, the recrossing correction of the Kramers-Langer approach becomes a large 
overcorrection and the resulting rate estimate too small by orders of magnitude. Analogous 
results were obtained by Shin et al. using one-dimensional Kramers’ theory in the globular 
limit pU] . 

The application of HTST to estimate the rate of escape of self-avoiding polymers is more 
challenging than the ideal polymer escape studied here. The reason is that the initial state 
where the polymer is coiled up corresponds to multiple local minima on the energy surface. 
The estimate of the initial state free energy, therefore, needs to take into account multiple 
conhgurations and configurational entropy rather than just the energy at a single minimum on 
the energy surface and vibrational entropy. One way to approach such multiple local minima 
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problems is systematic coarse graining |29j . Alternatively, the full TST can be applied as 
opposed to HTST where, for example, the reversible work formulation is used to identify an 
optimal transition state and estimate the free energy change in going from the initial state to 
the transition state [301 EH Eg. 
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Figure 1: The external potential of Eq. Q. The maximum of height AV = ~ 0.56 

is located at x = 0 and the minima are located at x = ±ao ~ ±1.22. The initial state, I, is 
conhned to the left well x < 0 and the right well x > 0 corresponds to the hnal state, F. In 
the multidimensional transition state for a polymer escape event, the center of mass of the 
polymer is located at the maximum of the external potential, as indicated by the label TS. 
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Figure 2: (a) The activation energy given by the maximum energy along the minimum en¬ 
ergy path on the energy surface as a function of N for iF = 10 (circles) and for K = 60 
(diamonds). The activation energy increases linearly with N until Nc, the onset of stretched 
polymer configurations, beyond which it remains constant, (b) Energy along MEPs vs. lo¬ 
cation of the center of mass of the polymer, Xq, for N = {1,16,64,120} K = 10 (circles) 
and for K = 60 (diamonds). The energy at the flat region becomes independent of Aq over 
a range of values for N > Nc 18 
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Figure 3: (a) The two lowest eigenvalues of the Hessian, Eq. (Iiup . as a function of chain 
length N. The negative eigenvalue A| corresponding to the unstable mode approaches zero 
as N increases beyond Nq- The smallest positive eigenvalue, A^, has a minimum of near 
zero at Nc and then reaches a constant value as N increases beyond Nc- (b) A saddle point 
conhguration and the modes corresponding to the two lowest eigenvalues of a polymer with 
= 56 beads and spring constant K = 60. The location of the beads is marked with black 
dots and the normalised eigenvectors indicated, si with triangle capped arrows and S 2 with 
square capped arrows. The si mode corresponds to motion along the minimum energy path 
towards the final state and the mode S 2 correfg)onds to the stretching of the polymer. 



































Figure 4: Rate of polymer escape as a function of chain length N for two values of the spring 
constant: (a) K = 10 and (b) K = 60. Circles and squares are results of direct simulations 
using Brownian and Langevin dynamics, respectively. Diamonds indicate rate estimates ob¬ 
tained from harmonic transition state theory with anharmonic corrections, HTST-I-AHC, for 
polymers with length close to Nq to remove divergence caused by the zero in the smallest 
positive eigenvalue of the Hessian at the saddle point. Stars indicate rate estimates obtained 
after correcting for recrossings, HTST-|-AHC-|-DC, using short time trajectories started at 
the transition state. The rate estimates obtained from hanger’s theory are indicated with 
triangles, showing a large underestimate of the rate for long chains. The eigenvalue A| is less 
than 10“® for N > 56 and K = 16 which causes problematic convergence of NEB along this 
mode leading to numerical inaccuracies in hanger’s rate. 
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Figure 5: The rate of escape of a polymer of length = 64 as a function of spring constant, K. 
Squares indicate the rates obtained from Langevin dynamics simulations, circles the HTST 
estimates, Eq. ([9]), stars the estimates from HTST with dynamical corrections, HTST+DC 
from (HSl), triangles the estimates from hanger’s theory, Eq. (dlD. The HTST+DC rate 
estimate is in good agreement with the rate obtained from the dynamics simulations, but 
hanger’s theory signihcantly underestimates the rate because the eigenvalue corresponding to 
the unstable mode approaches zero for long polymers. 
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Figure 6: The escape rate calculated for polymers with spring constant ol K = 60 as a 
function of length close to the crossover, Nc- Circles and squares indicate results obtained 
from dynamical simulations using Brownian and Langevin dynamics, respectively. Triangles 
indicate the HTST rate estimates, Eq. ([9]), diamonds the estimates from HTST with anhar- 
monic corrections, Eq. m, and stars the estimate with dynamical corrections, Eq. m- 
The anharmonic correction factor, g{a) given by Eq. (jl2p . shown in the inset, significantly 
reduces the diverging peak in the HTST rate estimate. 
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